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ABSTRACT 

The spatial distribution of matter in clusters of galaxies is mainly determined by the 
dominant dark matter component, however, physical processes involving baryonic mat- 
ter are able to modify it significantly. We analyse a set of 500 pc resolution cosmological 
<**j . simulations of a cluster of galaxies with mass comparable to Virgo, performed with the 

Q_i AMR code RAMSES. We compare the mass density profiles of the dark, stellar and 

gaseous matter components of the cluster that result from different assumptions for 
the subgrid baryonic physics and galaxy formation processes. First, the prediction of 
a gravity only N-body simulation is compared to that of a hydrodynamical simulation 
with standard galaxy formation recipes, then all results are compared to a hydrody- 
namical simulation which includes thermal AGN feedback from Super Massive Black 
Holes (SMBH). We find the usual effects of overcooling and adiabatic contraction in 
the run with standard galaxy formation physics, but very different results are found 
when implementing SMBHs and AGN feedback. Star formation is strongly quenched, 
producing lower stellar densities throughout the cluster, and much less cold gas is 
available for star formation at low redshifts. At redshift z = we find a flat density 
core of radius 10 kpc in both of the dark and stellar matter density profiles. We specu- 
late on the possible formation mechanisms able to produce such cores and we conclude 
that they can be produced through the coupling of different processes: (I) dynamical 
friction from the decay of black hole orbits during galaxy mergers; (II) AGN driven 
gas outflows producing fluctuations of the gravitational potential causing the removal 
of collisionless matter from the central region of the cluster; (III) adiabatic expansion 
in response to the slow expulsion of gas from the central region of the cluster during 
the quiescent mode of AGN activity. 

Key words: black hole physics - cosmology: theory - cosmology: large-scale structure 
of Universe - galaxies: formation - galaxies: clusters: general - methods: numerical 



1 INTRODUCTION 

Clusters of galaxies are the most massive virialised struc- 
tures observed in the Universe and provide a wonderful lab- 
oratory to test astrophysical theories. In the ACDM cos- 
mological scenario, clusters are assembled via a hierarchy of 
mergers of less massive structures like galaxies and groups of 
galaxies. Many physical processes play a role during the for- 
mation of a cluster. When satellite galaxies are accreted into 
a cluster, their properties can be changed by tidal and ram 
pressure stripping, leading to the formation of a wide variety 
of galaxy morphologies. Furthermore, clusters are known to 
be dark matter dominated structures, with most of the bary- 
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onic matter residing in a hot diffuse X-ray emitting gaseous 
phase, the Intracluster Medium (ICM). The stellar mass is 
less significant and mainly contained in the massive central 
elliptical galaxy. 

Since they are dominated by dark matter, this mass 
component determines the global properties of the mass dis- 
tribution in the cluster. However, from the theoretical side, it 
is well known that baryonic processes can produce significant 
differences in the distribution of matter in collapsed struc- 
tures with respect to models including only collisionless cold 
dark matter. For example, baryons are known to condense 
the centre of dark matter halos due to dissipative processes, 
produ cing adiabatic con traction of the total mass distribu- 
tion jGnedin et alj|2004h . Several models including baryonic 
physics have been invoked to solve the so-called cusp/core 
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problem in the ACDM cosmological framework, i.e. the dis- 
crepancy between the centrally cuspy dark matter profiles 
observed in dark matter halos in numerical N-body simula- 
tions and the centrally cored dark matter profiles inferred 
by observations in dwarf galaxies a nd low surface brightness 
galaxies (see the recent review bv lde Biokl linnl and refer- 
ences therein). The study of baryon physics induced modifi- 
cations in the mass distribution in collapsed structures, and 
clusters of galaxies in particular, is still a field with many 
open issues. 

The present paper is dedicated to the study of the ef- 
fects of baryonic processes on the mass distribution in clus- 
ters of galaxies. In particular, we use a set of cosmologi- 
cal hydrodyna mical simulation s performed using the AMR 
code RAMSES (|TevssierH200^ ) to study the effect of differ- 
ent models for baryons and galaxy formation physics on the 
mass density profile of a cluster of galaxies comparable to the 
Virgo cluster. This work can be considered a s an e xtension 
of the analysis performed by iTevssier" et all |201ll ) on the 
same sim ulations, and is compl ementary to the analysis pre- 
sented bv lMartizzi et all (2012). Here, we focus on the pecu- 
liar properties produced in the mass density profile when in- 
cluding Super Massive Black Holes (SMBHs) and the related 
AGN feedback in the recipes for galaxy formation physics. 
We stress that AGN feedback was initially introduced to 
solve the so-called "overcooling problem", namely the fact 
that too much stellar mass is produced in massive collapsed 
structure in hydrodynamical simulati ons with respect to 
what is observed in the real Universe (|Borgani fc Kravtsovl 
2009). The strong quenching of star formation produced by 
processes that couple AGN activity with the gas is expected 
to improve the match between simulations and observations 
llTabor fc Binnevll 19931 ; ICiotti fe Ostrikerll 19971 : ISilk fc Reesi 
1 199ST ) . Strong evidence for the existence of AGN feedback is 
provided by observations of X-ray cavities and radio blobs 
in galaxy clusters, typically interpreted as buoyantly rising 
bubbles of high entropy material injected in the central re- 
gion of clusters by jets of relativistic particles. In this paper, 
we show that by including SMBH physics and AGN feed- 
back it is possible to obtain interesting predictions on the 
modifications they can induce on the mass distribution in 
massive dark matter halos. 

The paper is organized as follows: the first section is 
dedicated to the numerical methods and the galaxy forma- 
tion recipes adopted for our simulations; we show our main 
results and provide our interpretation in the second section; 
the last section is left for a short summary of our results and 
to a discussion. 



2 NUMERICAL METHODS 

In this section, we describe the numerical techniques used 
to model the cluster we consider in this paper. We consider 
two cosmological hy drod ynamical simulations presented in 
ITevssier et all (|201lh and lMartizzi et all {{201$ ), plus a third 
cosmological simulation with only dark mat ter. For all th e 
three runs we used the AMR code RAMSES (|Tevssierll2002f ). 
The simulations were performed using the zoom-in technique 
which allows to obtain the required effective resolution only 
in selected regions of the computational domain. For all the 
three simulations the computational domain is a cubic box 



Cosmological parameters 



Type 


H [km s^Mpc" 1 ] 


CT8 


Q A 




DMO 


70 


0.77 


0.7 


0.3 


HYDRO 


70 


0.77 


0.7 


0.3 0.045 



Table 1. Cosmological parameters adopted in our simulations. 
The DMO label refers to the dark matter only run. The HYDRO 
label refers to the hydrodynamical runs. 



of side 100 Mpc/h. For the dark matter only run (dubbed 
as DMO, Table Q] and we adopt a standard ACDM cos- 
mology with parameters f2 m = 0.3, S1a = 0.7, ttb = 0.0, 
a s = 0.77 and H = 70 km/s/Mpc. For the two hydrody- 
namical runs (HYDRO runs in Table Q] and [5| we adopt the 
same values for f2 m , £Ia, 0"8 and Ho, but we set fib = 0.045. 
The initial condi tions for the thre e simu lations where com- 
puted using the lEisenstein fc Hul (I1998T) tr ansfer function 
and the grafic package 1 Bertschingerl 200ll) in its parallel 
implementation mpgrafic ( Prunet et alj 120*081 ). To perform 
the zoom-in technique we adopted the following approach: 
first, we ran a low resolution dark matter only simulation, 
then we identified dark matter halos at z = 0. From this 
halo catalog we built a set of halos whose virial masses lie in 
the range 10 14 to 2 x 10 14 M Q /h. Finally, we identified the 
final halo based on its assembly history: the halo is already 
in place at z = 1, therefore it can be considered as relaxed 
at z = 0. In particular the last major merger is observed at 
redshift z ~ 1.3. The final virial mass is M vir ~ 10 14 M , 
while M 20 oc = 1.04 x 10 14 M Q or M 50 oc = 7.80 x 10 13 M 0) 
where indice c refers to the critical density. This halo has 
been re-simulated at higher resolution, focusing the compu- 
tational resources in the region of the computational domain 
where it forms. 

In the DMO simulations the dark matter particle in the 
high resolution region has a mass of 9.6 x 10 6 M0/I1. In the 
two hydrodynamical runs, the dark matter particle in the 
high resolution region has a mass of 8.2 x 10 6 M /h, whereas 
the baryon resolution element has a mass of 1.4 xlO 6 M /h. 
Hydrodynamics is solved on an AMR grid that was initially 
refined to the same level of refinement than the particle grid 
(2048 3 , level t = 11). During the runs, 7 more levels of re- 
finement were considered (level ^ max = 18). The refinement 
criterion we used allows spatial resolution to be roughly con- 
stant in physical units; the minimum cell physical size was 
always close to Ai m i„ = L/2 £max ~ 500 pc/h. The grid 
was dynamically refined using a quasi-Lagrangian strategy: 
when the dark matter or baryons mass in a cell reaches 8 
times the initial mass resolution, it is split into 8 children 
cells. 

In the HYDRO runs, gas dynamics i s modeled us- 



ing a second-order u n split Godunov scheme l|Tevssierll2002l : 
ITevssier etafl 120061 : iFromang et al.l 120061 ) based on the 
HLLC Riemann solver and the MinMod slope limiter 
l|Toro et al.lll994l ). We assume a perfect gas equation of state 
(EOS) with 7 = 5/3. The HYDRO runs include standard 
subgrid models fo r gas cooling (accounting fo r H, He and 
metals; we use the ISutherland fc Dopitalll993l cooling func- 
tion), star formation (we adopt a star formation efficiency 
e» = 0.01) and super novae feedback (we adopt the "de- 
layed cooling" scheme, IStinson et al.l 120061 ) and metal en- 
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Mass and spatial resolution 



Type 


^cdm 










[10 6 M Q /h] 


[10 6 M /h] 


[10 6 M Q /h] 


[kpc/h] 


DMO 


9.6 


n.a. 


n.a. 


0.38 


HYDRO 


8.2 


1.4 


0.3 


0.38 



Table 2. Mass resolution for dark matter particles, gas cells and 
star particles, and spatial resolution (in physical units) for our 2 
sets of simulations. 



richment. In one of the two hydrodynamical simulations we 
als o implement AGN feedb ack, using a modified version of 
the lBooth fc Schavd (|2009l ) model. SMBHs are modeled as 
sink p articles, following the prescriptions of iKrumholz et all 
(2004), and they are allowed to merge. Gas accretion onto 
each SM BHs is computed adop ting a modified Bondi-Hoyle 
formula l|Booth fc Schavell2009l ). A fraction of the accreted 
mass is converted into thermal energy that is directly in- 
jected into the gas surrounding the black hole. For practi- 
cal purposes, in the rest of the paper we will refer to the 
run with AGN feedback as AGN-ON, and to the run with- 
out AGN feedback as AGN-OFF. Further details about the 
ga laxy formation and A GN feedback recip e s can be found 
in lMartizzi et al.l (|2012l) andlTevssier et~afl (|201lh . 

As discussed in iTevssier et all l|201lT ), the AGN feed- 
back energy is effectively deposited in two different modes: 
the energetic and impulsive quasar mode (Eddington lim- 
ited luminosity ~ 5 x 10 46 erg s" 1 ) during cold gas accre- 
tion when the accretion rate is high, and the quiescent radio 
mode during hot gas accretion when the accretion rate is 
low (Bondi-Hoyle- limited luminosity ~ 5 x 10 41 erg s _1 ). At 
z ~ the central SMBH in the AGN-ON cluster accretes 



mass in radio mode at a rate ~ 10 Mq yr , producing a 
total luminosity of 9 x 10 40 erg s _1 , comparable to the X-ray 
lumi nosity of the AGN in M87 Lx,o.5-7fc e v ~ 7 x 10 40 erg 
s _1 I Allen et al.l 120061 '); the mechanical power of the jet of 
the AGN in M87 is much larger, P w 3.3 x 10 43 erg s" 1 . 
These considerations imply that the AGN in M87 lies in an 
intermediate state between the radio and the quasar modes 
of our model. Further indications that the activity of the 
central SMBH is not particularly atypical are provided by 
comparison of the accretion rates and star formation rates in 
the AGN-ON simulation with observational results at red- 
shift z ~ 0. The mass accretion rate onto the SMBH is 
~ 10 -4 Mq yr -1 and the star formation rate within the 
inner few kpc from the cluster center is ~ 10 -2 Mq yr" 1 . 
T hese values are comparable with those recently measured 
bv lDiamond-Stanic fe Riekd l|2012h for Seyfert 1 galaxies. 



3 RESULTS 

In this Section we show how SMBHs and AGN feedback 
are able to change dramatically the properties of the mass 
distribution in clusters of galaxies. We focus on 3D mass 
density profiles of the cluster, analysing separately the dark 
matter and stellar components, as well as the total mass 
distribution. The goal is to highlight differences between the 
profiles in the DMO, AGN-ON and AGN-OFF run. The 



interpretation of our main results are provided in the next 
Section. 



3.1 Mass density profiles at redshift z = 

Spherically averaged mass density profiles of the cluster 
at redshift z = have been computed for all the sim- 
ulations. A reliable identification of halo centres is re- 
quired to compute density profiles. Halo centres and virial 
radii have been ide ntified using the AdaptaHOP algorithm 
llAubert et al.ll2004h. in the version imp emented and tested 
bv lTweed et al.1 (|2009r ). using the Most Massive Substruc- 
ture Method (MSM) to identify halos as well as their sub- 
structures. AdaptaHOP can be used to find the centres of 
groups of dark matter particles (i.e. "halos") as well as 
groups of star particles in the HYDRO runs (i.e. "galax- 
ies"). In the DMO run we directly use the centre of the dark 
matter halo, whereas in the AGN-ON and AGN-OFF runs 
we use the centre of the most massive group of star particles 
(i.e. the centre of the brightest cluster galaxy) since it traces 
better the minimum of the gravitational potential. To test 
that the latter choice does not lead to errors in the estimate 
of the profiles in the inner region of the cluster, we verified 
that picking the halo centre instead of the galaxy centre does 
not change the measured profiles. 

Figure[T]shows the dark matter and stellar mass density 
profiles at z = for the AGN-ON (solid line) and AGN-OFF 
(dashed line) runs; in this plot we use physical units. In the 
left panel the result of the DMO run (dotted line) is also 
plotted for comparison. All the profiles extend to the virial 
radius of the cluster. The three dark matter profiles look 
quite similar in the outer regions (r > 50 kpc), while they are 
significantly different in the inner regions. The AGN-OFF 
profile is much more centrally peaked and concentrated than 
the DMO profile, a result that can be interpreted as a result 
of adiabatic contraction of dark matter halos in response to 
the condensation of baryo ns at their centres l|Gnedin et al.l 
120041 ; I Tevssier et alj|2~01lf ); we test the quality of adiabatic 
contraction models in reproducing our results in the Subsec- 
tion !3.2l The AGN-ON profile is much different, with a very 
shallow inner slope. Within the inner 10 kpc we observe a 
dark matter core whose size is much larger than our spa- 
tial resolution limit (grey shaded area), whereas the AGN- 
OFF and DMO profiles seem to be consistent with cuspy 
dark matter profiles. Observational support to the results 
of our AGN-ON model comes from studies where gravita- 
tional lensing and dynamical data are combined to estimate 
the dark matter profiles of clus t ers of galaxies ijSand et al.l 
|2004 120081 ; iNewman et ai1l2009l , 1201 ll ; iRichtler et al.ll201ll l: 
these studies provide significative evidence of the existence 
of clusters of galaxies with dark matter profiles presenting 
cores or shallower inner slopes than in the Navarro-Frenk- 
White model. 

Looking at the stellar mass density profiles, we see that 
significantly less stellar mass forms in the AGN-ON run with 
respect to the AGN-OFF run that is affected by overcooling. 
At r > 50 kpc the stellar mass density in the AGN-OFF run 
is ~ 10 times bigger than in the AGN-ON run in all the ra- 
dial bins, and a very significative difference can be observed 
between the two profiles within ~ 10 kpc from the centre: 
the AGN-OFF density profile is increasing towards the cen- 
tre, whereas the AGN-ON profile shows a flat core with a 
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1 10 100 1000 1 10 100 1000 

r [kpc] r [kpc] 

Figure 1. Mass density profiles at z = 0. Left: dark matter. Right: stars. In all panels the grey shaded area represents the spatial 
resolution. 



size comparable to that observed in the dark matter profile. 
It is interesting that cores in the stellar surface brightness 
profiles of very massive elliptical galaxies and clu ster central 
galax i es have been observ e d by several authors (IKormendv 



19991 : iQuillen et ail l2000l; lEaine et all 120031: iTruiillo et al 



2004 ; Lauer et all 120051 ; ICote et alj|2007t iKormendv et al. 
2009l : lGrahamll201ll " and t hey can be deproiected in cored 
3D stellar density profiles l|Terzic fc Gr aham 20Q5|). 

In the AGN-OFF cluster the mass budget in the central 
region is strongly dominated by stellar mass: in absence of 
AGN feedback the star formation process is very efficient in 
turning gas into stars; the large stellar mass in the central 
region of the AGN-OFF cluster has been assembled through 
cooling flows providing gaseous material to trigger intense 
star formation events and through the accretion of satellite 
galaxies onto the central galaxy. In contrast, the mass in 
the central regions of the AGN-ON cluster has comparable 
contributions from stars and dark matter: overcooling is pre- 
vented from happening by AGN feedback and star formation 
is strongly quenched. 



3.2 Dark matter profiles at z = 0: testing the 
adiabatic contraction model 



It is possible to estimate some of the effects of baryonic pro- 
cesses, SMBHs and AGN feedback using an approach sim- 
ilar to that used by observers when analysing the surface 
brightness profiles of galaxies. First, an analytical model 
is used to fit the profile, then any significant deviation 
from the model is interpreted as a signature of physical 
processes. The same approach has been used to detect 
central light excesses/deficiencies with respect to a Sersic 
fit to the surface brightness profiles of early-type galax- 



ies (IKormendv et alll200alGrahamll201ll '). In lMartizzi et all 
|2012l ) we also adopted this approach to discuss the proper- 
ties of the stellar core observed in the stellar mass density 
profile in the AGN-ON run, showing that the Sersic function 
can be used to fit the stellar mass surface density profile out- 
side the cored region. Here, we use this criterion to analyse 
the dark matter profiles we measure in our simulations at 

2 = 0. 

We adopt the Einasto profile as our fiducial analytical 
model, since it has been shown to provide excellent fits to 
the dark ma tter profiles observed in cosmological N- body 
simulations (|Merritt et al. 21 11 1.",: iGraham et al. 2006). We 
use the following parameterization: 



jOEin(r) = p c exp < -d n 



where 



l/n 



d n = 3n - 1/3 - 0.0079/n. 



(1) 



(2) 



We use this analytical function to fit the dark matter profiles 
at z = 0, leaving I c , n and r c as free parameters. Our fits are 
performed using the Levenberg-Marqu ardt nonlinear least 
squares fit algorithm (|Press et al.ll 19921 ). The parameters we 
obtain after the fits are summarized in Table [3] In the DMO 
case we measure a concentration C200 = i?200c/-R-2 = 5.88, 
where R-2 is the radius at which the logarithmic slope of 
the dark matter pro file is -2. This val ue is typical for halos 
of mass ~ 10 14 M Q l|Reed et alll201lh . 

The fits are compared to the measured profiles in Figure 
[2] The Einasto profile provides a very good fit to the dark 
matter profile of our halo in the DMO run, despite the fact 
that it is typically used to fit the averge profile in cosmolog- 
ical simulations whereas we analyse only one halo. Both in 
the AGN-ON and AGN-OFF we observe that the Einasto 
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Figure 2. Comparison between the dark matter profiles measured at redshift 2 = (black solid lines) and the Einasto fits (red dashed 
lines). The scatter due to Poisson noise in each radial bin is also represented as a dark grey shaded area. The DMO Einasto fit is 
adiabatically contracted to get the two blue dashed profiles. Top-left panel: DMO run. Top-right panel: AGN-OFF run. Bottom-left 
panel: AGN-ON run. Bottom-right panel: relative difference between the measured profile and the Einasto fit as a function of radius (red 
lines); we use a dotted line for the DMO run, a dashed line for the AGN-OFF run and a solid line for the AGN-ON run; the blue lines 
show the residuals of the adiabatic contraction models with respect to the Einasto fits, fn all panels the grey shaded area represents the 
spatial resolution. 



profile is a good fit for r > 10 kpc, while we observe signi- 
ficative deviations with respect to the fitting formula in the 
inner regions (see the bottom- right panel of Figure [2}. The 
presence of these features can be interpreted as a manifesta- 
tion of processes that influence the formation of a standard 
distribution of dark matter in phase space; since these fea- 



tures are not observed in the gravity-only DMO run, we try 
to give them an explaination in terms of baryon physics. 

Dark matter halos are expected to respond adiabat- 
ically to the condensation of baryons at their centres 
( Gned in et al ] |2004l ; iTevssier et al.ll201ll ). We compare the 
results of our hydrodynamical runs at z = with the pre- 
diction of a simple adiabatic contraction model used by 



6 D. Martizzi et al. 



Einasto Profile Fits - z = 



Simulation 




Pe 


[M /fcpc 3 ] 


r e [kpc] 




n 


X 2 


DMO 


1.43 


X 


10 3 


± 1.5 x 10 2 


1.60 x 10 3 ± 6 x 


10 1 


5.93 ± 1.2 x 10- 1 


0.94 


AGN-ON 


1.39 


X 


10 3 


± 1.1 x 10 2 


1.55 x 10 3 ± 5 x 


10 1 


5.65 ±9 x 10" 2 


0.97 


AGN-OFF 


6.2 


X 


10 2 


±9 x 10 1 


2.0 x 10 3 ± 1.1 x 


10 2 


7.38 ± 1.8 x 10- 1 


2.21 



Table 3. Best fit parameters p c , r c and n for the Einasto profile at redshift z = 0. The value of the reduced chi-squared x 2 (1386 d.o.f.) 
for 10 kpc < r < 1000 kpc is also reported. The values are reported for the DMO, AGN-ON and AGN-OFF simulations. 



Adiabatic Contraction Model Parameters 



Simulation 


m d [M e ] 


r d [kpc] 


AGN-ON 


1.7xl0 13 


700 


AGN-OFF 


2.6X10 11 


10 



Table 4. Parameters adopted for the adiabatic contraction 
model: truncated sphere mass and radius . 

Abadi et al.1 (|2009^ . and already adopted in l|Tevssier et al.l 
20111 ); details for this model can be found in Appendix [Al 
We adiabatically contract the Einasto fit to the DMO pro- 
file assuming that the baryons are distributed in a constant 
surface density, truncated sphere of mass rrid and radius rd, 
despite the fact that the actual distribution is different. The 
values of nid and rd adopted in this paper are shown in 
Table I?] The result is plotted in Figure [2] (blue lines): in 
both the AGN-ON and AGN-OFF cases the adiabatically 
contracted profiles match the Einasto fits quite well at all 
radii, but, like the fits, the model does not provide a satisfac- 
tory description of the measured profiles for r < 10 kpc. A 
mass excess with respect to the adiabatic contracted profile 
is detected in the AGN-OFF case. The adiabatic contrac- 
tion model works reasonably well, especially at radii beyond 
a few kpc. In the interval 10 kpc < r < 1000 kpc we mea- 
sure reduced chi-squared values for the adiabatic contraction 
models y 2 = 2.14 for the AGN-OFF case and x" 2 = 0.72 for 
the AGN-ON case. These values are quite similar to what we 
obtain for the Einasto fits in the same radial range (Table 

EJi. 

In the AGN-ON run we find a mass deficiency with re- 
spect to the adiabatically contracted profile within 10 kpc 
from the centre; this mass deficiency is Mdcf = 5.7 X 10 10 
Mq . The fact that adiabatic contraction works for the outer 
regions of the cluster, but cannot predict the formation of 
the core, highlights the fact that additional processes play 
a role in shaping the properties of the mass distribution. 
Given the theoretical and computational limits of our phe- 
nomenological model for AGN feedback and the formation 
of SMBHs, and given that we limit our analysis to the evolu- 
tion of one halo, these results provide evidence that SMBH 
physics is indeed relevant in shaping the properties of clus- 
ters of galaxies, especially in the most overdense regions 
where SMBHs are expected to be found. We will discuss 
these topics in further detail in subection 13.41 

3.3 Evolution of the mass density profiles 

The study of the evolution of the mass density profiles in 
our simulations can be used to shed more light on the mech- 




1 10 100 1000 

r [kpc] 

Figure 3. Evolution of the dark matter density profile in the 
DMO run from z = 5 to z = 0. The grey shaded area represents 
the spatial resolution. 



anisms through which SMBHs and AGN feedback influence 
the mass distribution in clusters. To do so, we found the 
centres of the most massive progenitors of the cluster in 
each simulation, computed density profiles at redshifts z — 
1, 2, 3, 4, 5 and then compared to the result at z — 0. All 
the profiles shown in this subsection are plotted in physical 
units. 

First of all, we analyse the dark matter profile evolution 
in the DMO simulation (Figure |3}. In the outskirts of the 
cluster, the profiles change their slope; this transition marks 
the virial radius of the cluster. At redshift z = 5 the mass 
distribution in the halo is still quite different than at later 
times, with a shallower central slope and a radius ifeooc ~ 60 
kpc. At later times the inner part of the profiles (r < 20 kpc) 
reaches stability and evolution can be observed only in the 
outskirts where additional mass collapses into the halo, in 
fact the high r tail of the profile is constantly extending 
towards bigger distances from the centre. At redshift z — 
we measure that ifeooc ~ 1 Mpc. This kind of evolution is a 
clear exemple of stable clustering. 

In FigureUwe show the evolution of the dark and stellar 
mass density profiles in the AGN-OFF run. The effect of 
adiabatic contraction due to the condensation of baryonic 
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Figure 4. Evolution of the density profiles in the AGN-OFF run from z = 5 to z = 0. Left: dark matter. Right: stars. In all panels the 
grey shaded area represents the spatial resolution. 



mass at the centre of the cluster is already evident at redshift 
z = 5: the density within the inner 10 kpc is almost an order 
of magnitude larger than in the DMO case. Similarly to the 
DMO the dark matter halo is assembled, the dark 

matter profile in the inner region appears to be stable from 
z — 4 to z = 0, while additional dark mass appears to 
be accreted at higher radii. The stellar mass profiles evolves 
mantaining its shape, but increasing in amplitude. New stars 
are constantly formed and a significant increase of the stellar 
mass can be observed between all the considered snapshots. 



Finally, Figure [S] shows the same profiles for the AGN- 
ON run, highlighting some of the most interesting peculiari- 
ties of our model. The evolution of the dark matter profile in 
the outer regions of the cluster (r > 20 kpc) is very similar 
to that observed in the DMO case, but the inner region is 
strikingly different. The effects of the processes leading to 
the formation to the core observed at z = can be directly 
observed. At redshift z = 5 the dark matter profile appears 
to be contracted with respect to the DMO case, but much 
less than in the AGN-OFF case. Instead of mantaining the 
original slope, like in the AGN-OFF case, the inner part of 
the dark matter profiles becomes shallower with time, until 
a core is formed. A very similar behaviour is observed for 
the stellar mass density profile. While stellar mass is always 
increasing with time at radii r > 20, the inner profile evolves 
in a very peculiar way: at redshift z = 5 the stellar profile 
is very cuspy, but it evolves becoming shallower and shal- 
lower, despite the total stellar mass within 10 kpc from the 
centre increases. Between z — 1 and z = the stellar profile 
becomes almost completely flat for r < 10 kpc. 



3.4 Formation of central cores in the density 
profile 

The most evident peculiarity of the mass distribution in our 
AGN-ON simulation is the central core observed in the stel- 
lar and dark matter density profiles. From the plots (Figure 
[5} in the previous subsection we see that the dark matter 
core forms more gradually than the stellar core. At z > 1 
we see that the dark matter profile gradually evolves to- 
wards a cored configuration, while the stellar profile builds 
up mantaining its cuspy inner shape. The most interesting 
transition towards the formation of the two cores happens 
between redshifts z = 1 and z = 0: in this interval the dark 
matter profile becomes extremely fiat, while the cusp at the 
centre of the stellar profile is completely erased. The fact 
that this behaviour is not observed in the AGN-OFF simula- 
tion suggests that the cores in the AGN-ON dark and stellar 
matter profiles do not form naturally, but they are generated 
through external processes influencing the dynamics of col- 
lisionless matter. The aim of this subsection is to elucidate 
the mechanisms that lead to core formation in our simulated 
cluster. 

Several models have been proposed to produce cores 
in the density profile of a distribution of collisionless mat- 
ter starting from a cuspy profile; these models were origi- 
nally developed to explain the dark mat ter cores observed 
in gas rich dwarf galaxies (|de Blokll2010l ). but may be ap- 
plie d, in principle, to an y distribution of collisionless mat- 
ter. |E1-Zant et al.l (|200lh showed that sufficiently massive 
gas clumps can disrupt cusps through dynamical friction; 
this is not observed in our AGN-ON run since gas clouds 
are easily disrupted by AGN feedback. |E1-Zant et al.l (2004) 
showed that galaxies moving within the dark matter back- 
ground and transferring their orbital energy to the dark 
matter via dynamical friction may contribute to the forma- 
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Figure 5. Evolution of the density profiles in the AGN-ON run 
grey shaded area represents the spatial resolution. 
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2 = 5 to z = 0. Left: dark matter. Right: stars. In all panels the 



tion of cores. More recent ly, the very high resolution simula- 
tion of lNaab et all (|200Sh showed a similar effect: repeated 
minor dry mergers lead to decreases of the central stellar 
density concentration. This result has been confirmed by 
subsequent numerica l experiments, like those presented by 
lLaporte et al.l l|2012l 'l. The dense cores of accreted galaxies 
able to survive for a few crossing times modify the stellar 
mass distribution through dynamical friction. We emphasise 
that the same principle applies to the dark matter compo- 
nent, due to its collisionless nature. This effect is not ob- 
served in our AGN-OFF run because completely dry merg- 
ers are rare due to the high gas fractions in galaxies in this 
simulation. In the AGN-ON run, instead, a sign ificant frac- 
tion o f the gas is expelled from galaxies and the lNaab et al.l 
(2009) mechanism can be more efficient. Unfortunately, it 
seems to be challenging to produce an extended flat core 
like the one observed in our AGN-ON model only through 
repeated dry minor mergers. 

Alternative mechanisms to produce cores in dark mat- 
ter profiles involve purely gravitational processes related to 
SMBHs. In the context of a ACDM cosmology, where mas- 
sive structures form through the hierarchical mergers of less 
massive structures, the formation of SMBH binaries is an 
expected result. At the centre of collapsed objects SMBHs 
form binary pairs whose orbits decay as they transfer their 
orbital energy to collisionless matter via three-body interac- 
tions: the result is that collisionless matter can be expelled 
from the central regions via the gravitational slingshot effect 
and a core is for med. This process is usually r eferred to as 
SMBH scouring (|Milosavlievic fc Merritdl2003h and it is ex- 
pected to remove ~ 2 — 4 times the mass of the SMBH 
form ed after the binary completely decays (|Merritt et al.l 
120071 '). SMBH scouring is important on spatial scales from 
100 to 1 pc, thus it cannot be resolved in our simulations 
which has 1 kpc force softening. At the softening length the 
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Figure 6. Top: Evolution of the central SMBH observed at red- 
shift z = 0. Bottom: Evolution of the distance between the central 
SMBH observed at redshift z = and its closest neighbour. 



mass is completely dominated by dark matter and stars so 
that SMBH binaries cannot form. 

Despite the fact that SMBH scouring is not resolved 
in our simulations, they are able to resolve another process 
able to produce cores, namely SMBHs sinking to the very 
central region due to dynamical friction during mergers. The 
efficiency of this process has been extensively studied by 
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iGoerdt et all (|2010h and their results show that the orbital 
energy tranferred from the SMBHs to collisionless matter 
contributes to the formation of cores. The efficiency of this 
process is not as high as in the SMBH scouring, but it is able 
to produce mass deficiencies comparable to the mass of the 
infalling SMBHs. In our AGN-ON simulation the mass of 
the central SMBH at z = is M BH = 4.2 x 10 9 M Q , so this 
process by itself is not able to explain the mass deficiencies 
we observe in the central regions of the cluster: = 
5.7 x 10 10 M for dark matter and Mf e f = 3.04 x 10 10 
Mq for stars, measured with respect to an Einasto fit and a 
Sersic fit to the entire profile, respectively. 

Recently, the pure N-body simulations performed by 
iKulkarni fc Loebl |201ll ) showed that when multiple (two or 
more) SMBHs are present in a halo, the core formation pro- 
cess is much more efficient than during the inspiral of a single 
SMBH. Considering this enhanced core formation via dy- 
namical friction it is possible to create mass deficits of more 
than 5 times the total SMBH mass. The top panel of Figure 
[5] shows the mass growth of the central SMBH at the centre 
of the cluster at z = 0. The mass growth of this black hole is 
very smooth for z > 2, but it proceeds mainly through ma- 
jor mergers (sudden jumps in the black hole mass) at lower 
redshifts. The bottom panel of Figure [6] shows the distance 
between the central SMBH and its closest neighbour as a 
function of redshift, dubbed as dsn- Any sudden increase 
in the SMBH mass Mbh due to a black hole merger in the 
top panel plot can be detected as an increase in c/bh, be- 
cause after mergers the second neighbour becomes the first 
neighbour. The continuous variation of dBH tells us that at 
any redshift black holes are moving in the halo and losing 
orbital energy because of dynamical friction. A significative 
number of major mergers of very massive black holes hap- 
pens a t z < 1. This fact suggests that the lKulkarni fc Loebl 
(|2011T ) mechanism is particularly efficient at redshifts z < 1. 

SMBH infall is not the only process contributing to the 
formation of the core. In the AGN- ON simulation stron g 
AGN driven outflows are observed (|Martizzi et al.l I2012T ) . 
AGN feedback greatly increases the local temperature and 
entropy of gas, then the high entropy material is tran- 
ported out of the central region of the cluster through 
convective motions. These outflows modify the local grav- 
itational potential and may cause expansion of both the 
dark and stellar mass distribution. Similar processes have 
been observed in numerical simulations in which gas out- 
flows generated by supernovae feedback are used to pro- 
duce cores in the dark matter profiles of dwarf galaxies 
llNavarro et al ]| 19961; IGnedin fc Zhadl2002i;lRead fc Gilmorel 
120051 ; iGovernato et al] l201Ct IPontzen fc Governatd l201lF 7 
In the simulations performed by Navarro et alj ( 199rJ ). the 
mass outflows are simulated by growing and rapidly remov- 
ing an idealised potential from the centre of an equilibrium 
realisation of a dark matter halo, showing that the natural 
consequence is the formation of a core. The efficiency of this 
core formation process is oc M\[^ c R^J^ , where Mdi sc is the 
mass of the disc and i?disc is its scale radius. 

The calculations performed by IGnedin fc Zhad (|2002T ) 
show that this mechanism is extremely inefficient when 
a single supernova explosion is considered, however 
I Read fc Gilmorel l|2005l ) showed that repeated explosions fol- 
lowed by gas outflows and subsequent infalls are able to ac- 
count for the formation of dark matter cores. The recent nu- 




Figure 7. Variation of the gas mass enclosed in spheres of radius 
R = 2, 5, 10 and 20 physical kpc with respect to the total mass 
in the same region. The absolute value of the fluctuations can be 
larger than 1 because the variation of the gas mass is divided by 
the total mass after the outflow. 



mer ical experiments performed b y IGovernato et al] (|2010h 
and IPontzen fc Governatd |201ll ) confirm the efficiency of 
this mechanism in a fully cosmo logical context. In particu- 
lar, IPontzen fc Governatd (|201 ll ) suggest that gravitational 
potential fluctuations induced by supernovae driven outflows 
happening on a timescale shorter than the dynamical time 
cause the expansion of the dark matter distribution and the 
formation of a core. For this mechanism to be effective it is 
required to have gas outflows able to remove a significant 
fraction of the mass enclosed in the region where the core 
forms. Since we have multiple epochs of AGN driven gas 
outflows, this mechanism is likely to be active also in our 
AGN-ON run. This is indeed the case: Figure [7] shows the 
gas mass fluctuations induced by AGN burst driven outflows 
can be high fractions of the total mass in the central regions. 
Regions close to the central SMBH (radius r < R — 2 — 5 
kpc) present extreme mass fluctuations on short timescales. 
Within r < R = 10 kpc fluctuations become smaller, but can 
be as large as 10% of the enclosed total mass. At higher dis- 
tances from the centre the effect of these mass fluctuations 
is almost undetected. This means that potential fluctuations 
will be particularly strong only within the inner 10 kpc from 
the centre, that is the region where the core is observed. It 
is interesting to note that the amplitude of the mass fluc- 
tuations is typically higher at z > 1, however quite strong 
fluctuations are also observed at low redshift. 

Before drawing our conclusions, we carefully analyse 
what happens in the central region of the cluster. Figure 
[8] shows how the mass distribution in the central region of 
the cluster is influenced by the core formation mechanisms. 
We plot the evolution of the mass enclosed within 10 phys- 
ical kpc for all the different components: dark matter (solid 
lines), stars (dotted lines) and gas (dashed lines). The core 
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Figure 8. Time evolution of the total dark, stellar and gas mass 
enclosed within 10 physical kpc in the AGN-ON run. 



formation processes produce a decrease in the dark mass 
within 10 kpc from the centre from ~3x 10 11 M© at z = 5 
to ~ 2 x 10 11 Mq at z = 0. The stellar mass in the centre in- 
creases with time only down to z — 1, staying approximately 
constant until z = 0. At high redshift the star formation rate 
in the central region is high and the concentration of the 
stellar mass distribution is boosted by star formation events 
(Figure [9}; at redshift z < 1 the star formation in the cen- 
tral galaxy is strongly quenched by AGN feedback (Figure 
[9]), so the concentration of the stellar distribution cannot 
be boosted by strong star formation events, thus letting the 
core formation processes be effective. At z < 1 when AGN 
feedback becomes very strong because of the presence of very 
massive black holes, a slow decrease in the gas mass in the 
centre is observed: a gas mass ~ 10 11 M© is slowly removed 
from the central region before z = 0. This slow decrease in 
gas mass is expected to produce an adiabatic expansion of 
the total mass distribution, which will also contribute to the 
formation of a central core. Furthermore, the cooling time of 
hot gas within the inner 30 kpc of the cluster centre is ~ 1 
Gyr, suggesting that in its quiet mode the AGN can slowly 
eject the gas that rains down onto the centre from the inner 
cooling flow. 

The general picture we find from our analysis seems to 
show that the interplay between dynamical processes con- 
nected to SMBHs and the effects of AGN feedback on the 
star formation and the gas spatial distribution may provide 
an explaination to what is observed in our simulations. The 
dark matter core starts forming at redshift z = 4 — 5 due to 
AGN burst driven gas mass fluctuations on short timescales. 
The process goes on with very high efficiency down to red- 
shift z = 1. At redshift z < 1 the infall of very massive 
SMBHs contributing to the core formation process compen- 
sates for the smaller amplitude of the gas mass fluctuations. 
The final result is a flat dark matter core at redshift z = 0. 
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Figure 9. Average star formation rate in radial bins in the AGN- 
ON run. The grey shaded area represents the spatial resolution. 

The formation of the stellar core proceeds quite differently, 
with star formation compensating for the ejection of stellar 
material at z > 1, preventing the early formation of a core. 
At z < 1 the decrease in the star formation rate caused 
by AGN feedback finally allows the stellar core to form by 
redshift 2 = 0. 

Given our conclusions, we must stress that a clearer 
picture on the core formation problem has still to be drawn. 
The problem needs to be studied within a larger class of ha- 
los with different merger histories and masses, and idealised 
simulations are probably required to shed more light on the 
core formation process in clusters. 



4 SUMMARY AND CONCLUSIONS 

We used the results of a set of cosmological simulations to 
study the effect of baryons in massive dark matter halos, fo- 
cusing on phenomena involving SMBHs and AGN feedback. 
We use the zoom-in technique to simulate the formation of 
a galaxy cluster of mass comparable to Virgo at high spatial 
and mass resolution. In this work we analyse three different 
simulations: the first is a gravity only run with no baryons 
(DMO), the second is a hydrodynamical run with standard 
galaxy formation physics (AGN-OFF), the third is a hy- 
drodynamical run with standard galaxy formation physics 
plus SMBHs and AG N feedeback (AGN-ON) . We adopt a 
modified version of the lBooth fc Schavel l|2009T ) model to im- 
plement thermal AGN feedback in the AGN-ON run. This 
model allows to reproduce the self-regulated growth of black 
holes throughout the cosmic ages and it is tuned to repro- 
duce the observed Mbh — o relation. 

Our analysis focuses on the evolution of the 3D mass 
density profiles of the cluster from redshift z = 5 to redshift 
z = 0. We find a number of interesting differences between 
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our models that highlight the importance of accounting for 
all the physical processes taking place during the formation 
and evolution of the most massive bound structures in the 
Universe. Our results can be summarized in the following 
points: 

• At redshift z — the dark matter density profiles of the 
three simulations are consistent with each other at distances 
r > 10 kpc from the centre, but they differ significantly 
in the central region. The DMO profile is consistent with 
a standard Einasto model. The AGN-OFF profile shows a 
mass excess with respect to an Einasto fit due to adiabatic 
contraction of the baryons at the centre of the cluster. The 
AGN-ON profile presents a flat core of radius 10 kpc, i.e. 
a mass deficiency with respect to an Einasto fit. The dark 
matter core forms gradually from z = 5 to z = 0. We stress 
that dark matter cores in dark matter p rofiles have been 
claimed to be observed in s evera l clusters (ISand et al. | |2004 



2008!; iNewman et~ai1l2009l , l201ll ; iRichtler et all 



2011 



• The stellar density profile at 2 = is also very differ- 
ent between the AGN-OFF and AGN-ON runs. The AGN- 
OFF profile is very peaked at the centre, while the AGN- 
ON profile has a core of the same size as that of the dark 
one. Unlike the dark core, the stellar one forms between 
8 = 1 and z — 0. Note that, cores in the surface bright- 



ies have been observed by several authors (IKormendv 


199S 


Quillen et al.l 2000 


; Laine et all 20031: Truiillo et al. 


2004 


Lauer et alj 120051; 


Cote et al.l |200?I: Kormendv et al. 


2009 



Graham! 1201 IT ) , so our model is not in disagreement with 
observations and it predicts the existence of dark matter 
cores associated to stellar cores in massive ellipticals. Fur- 
thermore, in the AGN-OFF run the stellar density is larger 
than in the AGN-ON run at all radii. This is a clear ef- 
fect of overcooling of gas leading to enhanced (and too effi- 
cient) star formation. In the AGN-ON run, AGN feedback 
strongly quench es star formation and the overcooling prob - 



strongly quenches star formation and the overcooling prob- 
lem is avoided (|Tevssier et al.1 l201ll : iMartizzi et al.1 12OI2I) 



The result is also that less cool gas is available at the cen- 
tre of the AGN-ON cluster with respect to the AGN-OFF 
run: gas is heated by AGN feedback and carried away via 
convective motions and slow adiabatic expansion. 

• The core in the dark matter end stellar mass density 
profiles is a peculiar feature of our simulation that includes 
the physics associated with active galactic nuclei. In this pa- 
per, we suggest that the coupling of several mechanisms is 
responsible for their formation. First, SMBHs transfer part 
of their orbital energy to collisionle ss matter via dynam ical 
friction during dry galaxy mergers l|Goerdt et al.|[201(il ). es- 
pecially at redshift z < 1. Second, AGN driven gas outflows 
modify the gravitational potential in regions close to SMBHs 
with resulting ejection of collisionless matter from the cen- 
tral region of the cluster; subsequent gas outflows followed 
by the central 'revirialisat ion' of the central material a re ex- 
pected to produce cores (jPontzen fc GovernatdlioTll ): due 
to the stronger AGN feedback at z > 1, this mechanism is 
more effective at high redshift, but preserves an important 
role at low redshift. Third, the central hot gas slowly cools 
radiatively, falling onto the SMBHs in convective flows and 
is subsequently ejected impulsively; the slow loss of mass 
from the central region will produce an expansion of the 
inner mass distribution. 



Despite assertions that AGN feedback can affect cen- 
tral cluster dynamics, in fact our results seem to be com- 
plementary to the those obtained by other authors for 
less massive dark matter halos ([Governato et al.l 120101 : 



iPontzen fc Governatd l201ll : iMaccio' et al.1 l201lh . Related 



processes involving baryons seem to be active at the low 
and high mass end of the galaxy mass function, especially 
mechanisms involving gas outflows produced by feedback. 

Given the fact that we perform the analysis on zoom 
simulations of one cluster, we need to stress that our result 
needs support from a suite of dedicated cosmological simula- 
tion aimed at exploring the effect of SMBHs, AGN feedback 
and baryon physics in general on the most massive clusters 
in the universe. Furthermore, a number of numerical exper- 
iments is also needed to explore in detail the core formation 
processes that have been considered in our discussion, with 
particular attention on their coupling. The role of baryon 
m ass outflows in galaxy c l usters has been recently studied 
by iRagone-Figueroa et all |2012l ) using collisionless matter 
simulations where gas is modeled as a time varying exter- 
nal contribution to the gravitational potential; their results 
provide important support to ours. To conclude, we are con- 
vinced that our results are robust enough to assert that the 
implementation of AGN feedback and SMBHs in cosmolog- 
ical hydrodynamical simulations is an important ingredient 
for modeling massive clusters of galaxies. 
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APPENDIX A: ADIABATIC CONTRACTION 
MODEL 

The simplified mod el we adopt is based on that used in 
iTevssier et"ai1l|201ll} . If one defines the initial radius of each 
dark matter shell as r;, then an adiabatic contraction model 
is able to predict its value after contraction rf. In our case 
we adopt the transformation 

where Mi and M{ are the cumulative mass distributions be- 
fore and after contraction, respectively. The final cumulative 
mass distribution can be computed as 

M t = M dm (r f ) + M bar (r f ) = / dm Mi(n) + A/ bar (r f ). (A2) 

where MiirCj is the initial mass distribution in the DMO 
case, Mbar(7*f) is the baryonic mass distribution and Md m (Vf) 
is the adiabatically contracted dark matter distribution. The 
dark mass fraction is computed as /dm = 1 — ?rid /A/200 • Our 
aim is to recover the contracted dark matter profile Mdm(ff) 
given Mi(n) and Mbar(ff). We assume that the Mi(n) can 
be described by the fit to the Einasto profile we obtain for 
the DMO run at redshift z = 0: 

Mi(n) = 4tt f ' r 2 p E in(r)dr. (A3) 
Jo 
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The baryonic mass distribution is modeled as a constant 
surface density sphere with size rd and mass ma: 



We obtain the n value associated to each n solving numer- 
ically Equation I All and naturally obtain the adiabatically 
contracted dark matter profile Mdm(ff) using Equation IA2I 
Finally, we estimate the dark matter density profile after 
contraction as 




Pdm(r) 



1 dM dm (r) 



(A4) 



47rr 2 dr 



